###################################################################################
###################################################################################
####Appendix: Non-state Atrocities in Capital Cities - Negative Binomial Models####
###################################################################################
###################################################################################
library(MASS) 
library(pscl)
library(foreign)
library(stargazer)
library(mvtnorm)
library(plyr)

# Set working library
setwd("~/Data/Global Analysis/")
# Read in main data
main.data <- read.dta("full.grid.upd.dta")

######################
###Global NB Models###
######################
###Model A-NB1 (corresponding to Models 1 and 2)
nb1 <- glm.nb(incidentnonstatefull ~ lag_Capital + lagcivconflagtemp + lagincidentsumfull
              + loglagppp + loglagbdist1 + lagp_polity2 + loglagpop + loglagttime + loglagcellarea + year96 + year97+
                year98 + year99 + year00 + year01 + year02 + year03 + year04 + year05 + year06 + year07 + year08, 
              data = main.data)
summary(nb1)
AIC(nb1)

###Model A-NB2 (corresponding to Model 3)
nb3 <- glm.nb(incidentnonstatefull ~ lag_Capital + lagcivconflagtemp + lagincidentsumfull + lagurban
              + loglagppp + loglagbdist1 + lagp_polity2 + loglagpop +  loglagttime + loglagcellarea + year96 + year97+
                year98 + year99 + year00 + year01 + year02 + year03 + year04 + year05 + year06 + year07 + year08, 
              link = log, data = main.data)
summary(nb3)
AIC(nb3)

##Model A-NB3 (corresponding to Model 4)
nb4 <- glm.nb(incidentnonstatefull ~ lag_Capital + lagcivconflagtemp + lagincidentsumfull + lagurban + loglagross_oil_prod
              + loglagppp + loglagbdist1 + lagp_polity2 + loglagpop +  loglagttime + loglagcellarea + year96 + year97+
                year98 + year99 + year00 + year01 + year02 + year03 + year04 + year05 + year06 + year07 + year08, 
              link = log, data = main.data)
summary(nb4)
AIC(nb4)

####Export to LaTex
stargazer(nb1, nb3, nb4)
AIC(nb1, nb3, nb4)

############################
###Urban Sample NB Models###
############################

###Model A-NB4 (corresponding to Models 5A-5C)
nb5 <- glm.nb(incidentnonstatefull ~ lag_Capital + lagcivconflagtemp + lagincidentsumfull + lag.large.urb + loglagross_oil_prod
              + loglagppp + loglagbdist1 + lagp_polity2 + loglagpop + loglagttime + loglagcellarea + year96 + year97+
                year98 + year99 + year00 + year01 + year02 + year03 + year04 + year05 + year06 + year07 + year08, 
              data = u.main.data)
summary(nb5)
AIC(nb5)

###Model A-NB5 (corresponding to Models 6A-6C)
nb6 <- glm.nb(incidentnonstatefull ~ lag_Capital + lagcivconflagtemp + lagincidentsumfull + lag.large.urb + loglagross_oil_prod
              + loglagppp + loglagbdist1 + lagp_polity2 + loglagpop + loglagttime + loglagcellarea + year96 + year97+
                year98 + year99 + year00 + year01 + year02 + year03 + year04 + year05 + year06 + year07 + year08, 
              data = u.main.data.Af)
summary(nb6)
AIC(nb6)

#####Export to LaTex
stargazer(nb5, nb6)
AIC(nb5, nb6)


########See Global and Urban dofile for comparison with ZINB models and Vuong non-nested hypothesis test results
